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ABSTRACT 

We propose a model of hard X-ray flares in protostars observed by ASCA 
satellite. Assuming that the dipole magnetic field of the protostar threads the 
protostellar disk, we carried out 2.5-dimensional magnetohydrodynamic (MHD) 
simulations of the disk-star interaction. The closed magnetic loops connecting 
the central star and the disk are twisted by the rotation of the disk. As the 
twist accumulates, magnetic loops expand and finally approach to the open field 
configuration. A current sheet is formed inside the expanding loops. In the 
presence of resistivity, magnetic reconnection takes place in the current sheet. 
Outgoing magnetic island and post flare loops are formed as a result of the 
reconnection. The time scale of this 'flare' is the order of the rotation period 
of the disk. The released magnetic energy partly goes into the thermal energy 
and heats up the flaring plasma up to 10 8 K. The length of the flaring loop 
is several times of the radius of the central star, consistent with observations. 
The speed of the hot plasmoid ejected by the reconnection is 200 — 400 km s _1 
when the footpoint of the loop is at 0.03 AU from 1 M Q protostar. The hot 
plasma outflow can explain the speed and mass flow rate of optical jets. Dense, 
cold, magnetically accelerated wind (v ~ 150 — 250 km s^ 1 ) emanates from the 
surface of the disk along the partially open magnetic field lines threading the 
disk. This dense, cold wind may correspond to high velocity neutral winds. 

Subject headings: stars: activity — stars: pre-main sequence — ISM: jets and 
outflows — accretion, accretion disks — MHD 
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1. INTRODUCTION 

Observations by ASCA revealed that protostars are strong hard X-ray emitting sources. 
By using ASCA, Koyama et al. (1994,1996a,b) carried out a systematic survey of hard 
X-ray sources in molecular clouds. They detected hard X-rays (4-10keV) associated with 
Class I IR sources (protostars). A protostar (Rl in Coronet cluster) exhibited a powerful 
X-ray flare whose spectrum showed significantly higher temperature (T ~ 10 8 K) than 
X-ray flares in T Tauri stars which is usually ~ (1 — 2) x 10 7 K (Koyama et al. 1996b). The 
X-ray luminosity in the 2-10keV band is ~ 10 31 erg s _1 and the e-folding time of the flare 
is ~ 5 x 10 4 sec (Koyama et al. 1996b). Three X-ray emitting class I sources in Rho Oph 
cloud also showed significantly hard X-rays. One of them showed flare like time variabilities 
(Koyama et al. 1994). 

Although the analogy with solar flares may be appropriate for such X-ray flares in 
protostars, they differ from solar flares in the following points: (1) The total energy released 
by the protostellar flare (10 35 — 10 36 erg) is 10 5 times larger than solar flares. (2) The 
spectrum of the X-ray emission in protostellar flares is harder than solar flares, indicating 
the existence of 8 keV plasma, and (3) The size of the X-ray emitting region is the order 
of 10 11 cm (Koyama 1995 private communication). If the emission comes from magnetic 
loops, their size should be several times larger than the protostar. 

It has long been recognized that solar flares are often associated with sheared magnetic 
fields. Mikic & Linker (1994) and Linker & Mikic (1995) carried out MHD simulations of 
the dynamical evolution of magnetic arcades in spherical geometry under shearing motions 
imposed on the line-tied footpoints of an initially current free arcade. When the amount 
of the footpoint shear is sufficiently large, the magnetic arcades begin to inflate. In the 
presence of finite resistivity, magnetic nonequilibrium leads to the magnetic reconnection 
and plasmoid ejection in the current sheet developed inside the expanding arcades. 
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In the innermost region of protostellar disks, the magnetic fields of the protostar will 
affect the dynamical evolution of the disk (e.g., Camenzind 1990). Since the rotation 
of the disk continuously injects helicity into the closed magnetic fields connecting the 
central star and the disk, we expect magnetic nonequilibrium and resulting magnetic 
reconnection similar to those in solar flares (see also Ghosh & Lamb 1979, 1991; Zylstra 
1988; Lynden-Bell & Boily 1994). Recently, the magnetic interaction between a central star 
and its surrounding disk has attracted a renewed interest in relation to the activity of young 
stars (Shu et al. 1994, Hirose 1994, Lovelace et al. 1995). In order to explain why T Tauri 
stars have rotation rates well below the break-up despite the high accretion rate, Shu et al. 
(1994) developed a steady "X-wind model" which emanates from the point where the disk 
is corotating with its central star. Through the magnetic interaction between the central 
star and the inner edge of the disk where they assumed to be the corotation point, excess 
angular momentum brought in through the disk are extracted by X-winds. Subsequently, 
Ostriker & Shu (1995) analyzed a steady state configuration of the magnetosphere between 
the central star and the disk and discussed how dipole fields become open and generate the 
"helmet streamer" by magnetic reconnection. 

In the steady X-wind model, Shu et al. (1994) assumed that the dipole magnetic fields 
of the central star do not thread the disk except the inner edge where they assumed an 
enhanced magnetic diffusivity. It is possible, however, that the magnetic fields of the central 
star thread the disk either by enhanced turbulent magnetic diffusivity inside the disk or 
by the recoupling of the magnetic field and the disk matter due to the enhanced ionization 
rate near the central star. 

In this Letter, we present the results of nonlinear, time-dependent MHD simulations of 
the interaction between the dipole magnetic field of the central star and the rotating disk 
threaded by the dipole field. 
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2. NUMERICAL MODELS 

We solve the resistive MHD equations in cylindrical coordinate (r, 0, z) by applying 
a modified Lax-Wendroff scheme with artificial viscosity (see Shibata 1983 and references 
therein). Axial symmetry is assumed. The effects of viscosity, radiative cooling and rotation 
of the central star are neglected. We take the z-axis parallel to the rotation axis of the disk 
and assume that the magnetic dipole moment is aligned with the rotational axis. 

Exact solutions of the disk in rotational equilibrium are obtained under the simplifying 
assumptions for the distributions of angular momentum L = L r a and pressure P = Kp 1+1 ^ n 
where we adopted a = 0.494 and n — 3. The equation of force balance is integrated to give 



(r 2 + z z ) l ' z 2(1 — a) p r 2(1 — a) p 

where M is the mass of the central star and the subscript denotes quantities at the 
pressure maximum (r ,0,0). The mass distribution inside the disk is determined by 
specifying the parameters a, n and Po/Po- We use the normalizations r = vko — Po — 1 
where v K0 = (GM/tq) 1 / 2 is the Keplerian rotation speed at r = r . The unit of time is 
T~o =r o/ v Ko- The mass distribution outside the disk is assumed to be that of the spherical 
isothermal corona surrounding the central star. The corona is assumed to be static at t — 0. 
The coronal density at (0, 0, r ) is p h = 2 x lCT 5 p - The coronal temperature T h is 7^=1.667 
xl0 2 T , where T is the disk temperature at pressure maximum. The total simulation 
region in the r — z plane is taken to be < r < 13.1ro and < z < ll.Oro- The number of 
the mesh points used in our calculation is (N r ,N z ) = (305,364). The grid spacings increase 
with r and z. 

We use free boundary conditions at the outer boundaries. At the equatorial plane 
(z = 0), we assume a symmetric boundary condition. We introduce a thin shell-like 
damping region (stellar surface) between r in = (r 2 + z 2 ) 1 / 2 = 0.336r and r s = 0.528r 
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where all quantities except magnetic fields are damped to the initial values. Alfven waves 
can be reflected in this layer. All quantities are set to the initial values in stellar interior 
(f < r in)- Meanwhile magnetic fields are stress free at the outer boundaries. 

As for resistivity, we assume anomalous resistivity r)(r,z) = r) (vd/v c — 1.0) 2 when 
v d = \ J\/p exceeds the critical value v c where J is the current density. When vj < v c , r)(r, z) 
is set zero. This type of anomalous resistivity has been used in simulations of magnetic 
reconnection between the emerging flux and the overlying magnetic fields in the solar 
corona (see e.g., Shibata, Nozawa & Matsumoto 1992, Yokoyama & Shibata 1995). 

3. NUMERICAL RESULTS 

The results of the simulation of a typical model are shown in Figure 1 and Figure 
2. The model parameters are A 1 = v 2 s0 /(^v 2 m ) = 2.0 x 1(T 3 , A 2 = v 2 A0 /v 2 K0 = 2.0 x 10~ 4 
where v s o and vao are the sound speed and the Alfven speed at (r , 0,0), respectively. 
The adiabatic index 7 is 5/3. We use the anomalous resistivity with v c = 1000 and 
i] = 0.01 in non-dimensional units. The plasma f3 in the disk at the pressure maximum is 
Po = 8nP /B 2 = 2A ± /A 2 = 20. 

Solid curves in Figure 1 show the magnetic field lines projected onto the r — z plane. 
The color scale shows the temperature. The arrows show velocity vectors in the r — z 
plane. The magnetic field lines connecting the central star and the disk are twisted by 
the rotation of the disk. As the magnetic twist accumulates, the magnetic loops begin to 
inflate. They expand quasi-statically in the early stage but later (t > 2r ), they expand 
dynamically. A current sheet is formed inside the expanding loops and the current density 
exceeds the threshold for the onset of anomalous resistivity. Around t = 2Atq, magnetic 
reconnection takes place in the current sheet. Outgoing magnetic island and post flare loops 
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are formed after the reconnection (t = 2.68r ). Hot plasmas created by the Joule heating 
and compressional heating are confined in the outgoing magnetic island. The speed of 
outflowing hot plasma is the order of the local Alfven speed va ~ (2 — 5)vko- Crest-shaped 
fast MHD shock waves are formed in front of the expanding magnetic loops. By parametric 
survey, we confirmed that when the coronal density is lower the speed of reconnection flow 
and temperature are higher because local Alfven speed increases. 

Inside the magnetic island, magnetic field lines are almost toroidal. Since they are 
detached from the star and the disk, the mass and the total angular momentum of the 
rotating toroid is conserved. By equating the centrifugal force l 2 /r 3 where / is the specific 
angular momentum with the radial component of gravity GMr/(r 2 + z 2 ) 3 ^ 2 , we get r/z — > 
as z — > oo ; the plasmoid is collimated toward the rotational axis. The azimuthal magnetic 
fields also helps collimating the outflow by self-pinching effect (see Uchida & Shibata 1985; 
Shibata & Uchida 1986). Although the computational domain in our simulation is too small 
to show the collimation, the velocity vectors at t = 4.0to show a tendency to point toward 
the polar directions. 

The magnetic field lines near the disk have small angle with respect to the equatorial 
plane. This configuration is favorable for the appearance of magneto-centrifugal winds from 
accretion disks (e.g., Blandford & Payne 1982). The speed of the disk wind is the order 
of the Keplerian rotation speed at the footpoint of the field lines on the disk. Since the 
density is high, the mass outflow rate of the cold disk wind is several times larger than that 
of the hot plasma jet driven by the reconnection. 

Figure 2 shows the density distribution (gray scale) and isocontours of magnetic energy 
density. Just after the reconnection, the density of the reconnection region decreases due to 
the heating and expansion. During the later evolution, the density increases to p ~lCr 5 p 0) 
because the colder gas that flows out of the disk refills the region. This justifies rather high 
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initial coronal density (ph ~ 2 x 10 5 p ) in our model. 

Figure 3 shows the time history of various energies. The solid curve in the lower panel 
shows the total magnetic energy E mag in the simulation region. The dotted curve shows the 
magnetic energy E h;mag of hot plasma (T > 1.2T h ). The magnetic energy increases almost 
linearly with time until magnetic reconnection takes place. The kink around t = 2.4r is 
due to the magnetic reconnection. The magnetic energy of hot component saturates and 
slightly decreases. In the MHD simulations of solar coronal magnetic loops which subject to 
the helicity injection, Kusano, Suzuki & Nishikawa (1995) noted similar kink and explained 
it by extending the Woltjer & Taylor's minimum energy principle (Woltjer 1958, Taylor 
1974). 

The solid curve and the dotted curve in the top panel of Figure 3 show thermal energy 
Eh,th, and the kinetic energy E hjkin of hot plasma (T > 1.27^), respectively. Rotational 
energy is included in Eh t kin- The kinetic energy slowly increases during quasi- hydrostatic 
evolution (t < 2r ). The kinetic energy increases more steeply when the equilibrium is lost 
(t ~ 2.0t ). Since the magnetic loops expand rapidly, the rate of increase in the coronal 
thermal energy E hyth decreases during t — 2.0 — 2.4t . The coronal thermal energy increases 
more steeply after magnetic reconnection at t — 2Atq. These curves clearly show the 
bifurcation of the system into a different energy state after magnetic reconnection. 

It should be noted that the energy injected into the corona through the rotation of the 
disk (E in j = E mag + E h th + E h>kin ) increases almost linearly with time even after magnetic 
reconnection because magnetic helicity is continuously injected from the disk. At t — 3.0r , 
about 70% of the injected energy goes into the magnetic energy, 18% to the thermal energy, 
and 12 % to the kinetic energy. The thermal energy of hot plasma at t — 3.0r is about 
40% of the initial total magnetic energy (15% of the total magnetic energy at t — 3r ). 
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4. DISCUSSION 

In this Letter, we demonstrated that if the dipole magnetic field of the central star 
is connected to the rotating disk which is not co-rotating with the central star, the 
dipole magnetic fields can disrupt by imposed twists within one rotation of the disk. We 
schematically summarize our numerical results in Figure 4. Our 2.5D nonlinear MHD 
simulation confirmed the previous expectation from the equilibrium calculations by Zylstra 
(1988) that magnetic reconnection will be triggered by the loss of equilibrium due to the 
helicity injection. We also showed that hot plasmoids are ejected in bipolar directions 
with velocity 2—5 times the Keplerian rotation speed around the inner edge of the disk. 
Plasma heating occurs both by the Joule heating in the current sheet developed inside 
the inflating magnetic loops and by shock waves created by magnetic reconnection. The 
plasma temperature can rise up to the temperature determined by v s ~ v a ~ 2 — 5vk- 
Since v K ~ 170 km s -1 at r = 0.03 AU around 1 M star, the temperature of the plasmoid 
and the post flare loop can be 10 7 — 10 8 K. According to the observation of an X-ray flare 
from a protostar (Koyama et al. 1996b), the emission spectrum of protostellar flare extends 
up to lOkeV. This result is consistent with the above plasma temperature expected from 
magnetospheric flares in protostars. 

Observation of X-ray flares by ASCA (Koyama et al. 1996b) showed that the e-folding 
time of the flare is ~ 5 x 10 4 sec and the total energy of the flare is 10 36 erg. The size of 
X-ray emitting region derived from the emission measure (~ 10 54 cm~ 3 ) and flare decay 
time is the order of 10 11 cm. The number density of the flaring plasma is nfi are ~ 10 10 5 
cm -3 . This rather high density can be explained by the dynamical refillment of the corona 
by the mass outflow from the disk (see Figure 2). When the X-ray emitting plasmas are 
confined in loops, their size should be larger than the radius of the protostar. Our numerical 
results are consistent with this observational constraint because hot plasmas are created in 
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gigantic loops several times larger than the stellar radius. The mass flux of the hot plasma 
outflow obtained from simulation is 

V0.03AU; U70kms-V V3-10 10 cm-V &y 
After the X-ray emitting gas is cooled down (cooling time ~ 1 day), it will be observed as 
optical jets. The mass flux of the hot plasma is consistent with the observation of optical 
jets. 

Next, let us discuss whether the magnetospheric flares can release the observed total 
energy of X-ray flares. In our simulation, initial magnetic energy is thermalized within one 
rotation period (~ 10 5 sec at 0.03 AU). If the size of the flaring region is d ~ 10 11 cm and 
the mean strength of magnetic field in the flaring region is B ~ 100 — 200 G, the released 
energy L ~ (B 2 /8ir) d 3 ~ 10 36 erg is the order of the observed luminosity. 

We checked the dependence of numerical results on resistivity. When the threshold for 
anomalous resistivity v c is larger, the plasma temperature and the speed of reconnection 
flow are larger because magnetic reconnection does not take place until large amount of 
magnetic energy is accumulated. Since the speed of the jet is larger, the magnetic island 
is deformed into a "heart shape" as shown in Figure I. This feature has also been noted 
by Choe (1995) when he simulated the solar magnetic arcades under photospheric shearing 
motion. When the resistivity is uniform or v c is small, the magnetic island is more circular. 

We would like to comment on the relation between our model and the X-wind model by 
Shu et al. (1994). They considered that the magnetic fields of the central star can thread 
the disk only in the prescribed small diffusive transition zone near the point where the disk 
corotates with the central star. Interior to this radius, magnetic field lines are closed and 
connected to the surface of the star. Along the open magnetic field lines, X-winds emanate 
from the small diffusive zone around the inner edge of the disk. This configuration is similar 
to the final stage of our simulation in which the magnetic field lines after reconnection 
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become partially open, and dense, cold winds flow out from the disk along the partially 
open field lines. Although we cannot continue the simulation until steady state is reached, 
we could show some essential components of the X-wind model. Of course, we need to 
include the effect of the stellar rotation in order to study the angular momentum exchange 
between the star and the disk. We want to extend our study in such direction in near 
future. 
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visit to the National Astronomical Observatory, Japan. We also thank Prof. F. Lamb for 
pointing out the related work by Zylstra (1988). Numerical computations were carried out 
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Fig. 1. — The result of the numerical simulation for a model with A\ = ^/(t^Ito) = 
2.0 x 10- 3 ,A 2 = v 2 A0 /v 2 K0 = 2.0 x 10" 4 and p h /p = 2.0 x 10~ 5 . The color scale shows 
the temperature normalized by Tk=pv 2 KQ / R, here \i is mean molecular weight and R is gas 
constant. The initial coronal temperature is T h /T k —l/3. Solid curves denote magnetic field 
lines. The region of closely spaced field lines is included to identify the reconnection region 
in the subsequent frames. Arrows depict velocity vectors in r — z plane. The unit velocity 
is shown by arrows at the top right of each panel. 

Fig. 2. — Time variation of the density distribution. Color scale shows the density. Solid 
curves denote isocontours of magnetic energy density. Arrows depict the velocity vector. 

Fig. 3. — Time history of various energies normalized by the initial magnetic energy. Upper 
panal shows the thermal energy Eh,th (solid curve) and kinetic energy Eh,kin (dotted curve) 
of hot component, 

E h , th = 2n / [P/( 7 - l)]rdrrfz, E hMn = 2ir / (pv 2 /2)rdrdz. 
Jo Jr e Jo Jt s 

The lower panel shows the total magnetic energy E mag (solid curve) and that of the hot 
component Eh tmag (dotted curve). 

Fig. 4. — A schematic picture of numerical results. Hot plasma jet ejected from the flaring 
region corresponds to the optical jet. Cold, dense wind emanating from the disk may explain 
high velocity neutral winds. 
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